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Abstract 

Cellular electrophysiology is often modeled using the cable equations. 
The cable model can only be used when ionic concentration effects and 
three dimensional geometry effects are negligible. The Poisson model, in 
which the electrostatic potential satisfies the Poisson equation and the 
ionic concentrations satisfy the drift-diffusion equation, is a system of 
equations that can incorporate such effects. The Poisson model is unfor- 
tunately prohibitively expensive for numerical computation because of the 
presence of thin space charge layers at internal membrane boundaries. As 
a computationally efficient and biophysically natural alternative, we in- 
troduce the electroneutral model in which the Poisson equation is replaced 
by the electroneutrality condition and the presence of the space charge 
layer is incorporated in boundary conditions at the membrane interfaces. 
We use matched asymptotics and numerical computations to show that 
the electroneutral model provides an excellent approximation to the Pois- 
son model. Further asymptotic calculations illuminate the relationship of 
the electroneutral or Poisson models with the cable model, and reveal the 
presence of a hierarchy of electrophysiology models. 

1 Introduction 

Electrophysiology is the study of the electrical activity of biological tissue [TJ [J] . 
Because of its importance in many physiological processes and its quantitative 
nature, it has been a favorite subject in biophysics and mathematical physiology. 
Traditional mathematical models of cellular electrical activity are based on the 
famous work of Hodgkin and Huxley [5] , and may be collectively termed cable 
models. These models are based upon an ohmic current continuity relation on 
a branched one dimensional electrical cable [T2J QT] . The derivation of the cable 
model is based on several important assumptions |12| : 

• A one dimensional picture, or more generally, a one dimensional tree rep- 
resentation of cell geometry is adequate. Geometrical details that are lost 



in making this simplified description have negligible effect on electrophys- 
iology. 

• The extracellular space can be reduced to a single isopotential electrical 
compartment. 

• Ionic concentrations are effectively constant in space and time within each 
cell separately and in the extracellular space. The diffusive current that 
may be induced by concentration gradients or the changes in equilibrium 
potential arc negligible. 

Such assumptions are justified in many instances, for example in the isolated 
neuronal axon [5], where the cable model has been extremely successful in ex- 
plaining the physiology and in making quantitative predictions - a triumph 
counted among the greatest successes of mathematics in biology. There may, 
however, be many cases in which any or all of the above assumptions are vio- 
lated especially in the central nervous system and cardiac tissue, as suggested 
by the complex microhistological structure they exhibit [HI H] ■ One line of work 
that addressed this difficulty was that of Qian and Sejnowski [20]. Their work 
addresses the last of the above difficulties, but retains the one-dimensional char- 
acter of the cable model. 

In |15| , we presented a three-dimensional model of cellular electrical activity 
which addresses all of the above limitations of the cable model. This model 
consists of a system of partial differential equations to be satisfied by the ionic 
concentrations and the electrostatic potential. In this paper, we introduce a 
slight modification of this model, which we call the electroneutral model. 

The first goal of this paper is to demonstrate the validity of the electroneutral 
model by comparing this with the Poisson model [T3] • In the Poisson model, the 
ionic concentration dynamics is governed by the drift-diffusion equations and the 
electrostatic potential satisfies the Poisson equation. Non-dimcnsionalization re- 
veals the presence of multiple temporal scales and of a thin boundary layer at 
the membrane interfaces in which electric charge accumulates (Debye layer) [53] . 
This boundary layer necessitates the use of a fine spatiotemporal mesh in numer- 
ical simulations making such computations prohibitively expensive. We intro- 
duce the electroneutral model as an alternative to the Poisson model, in which 
the Poisson equation is replaced by the electroneutrality condition. The model 
does not resolve the dynamics within the thin boundary layers and instead in- 
corporates the effect of these layers by modifying the boundary conditions at 
the membranes. The boundary layers are incorporated as charge densities of 
zero thickness at the membrane, a picture that is better aligned with the bio- 
physical view of the membrane being a capacitor within a conducting medium. 
This obviates the necessity for high spatiotemporal resolution in computations, 
making the electroneutral model far more amenable to numerical simulation. 
Using matched asymptotics, we show that the electroneutral model provides 
an approximation to the Poisson model. We present computational studies in 
the final section to demonstrate that the electroneutral model does indeed pro- 
vide an excellent approximation to the Poisson Model for biophysically realistic 
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Figure 1: The variables <f>, Cj are defined in the regions 

an d QW, which 

we have identified as intracellular and extracellular regions in the above. The 
membrane acts primarily as a capacitor, but possesses ionic channels through 
which transmembrane current can flow. 



parameter values. 

The second goal of this paper is to clarify the relationship between the 
Poisson and electroneutral models to cable models. If we are to claim that the 
Poisson or electroneutral models arc a generalization of the cable model, we 
would like to know under what conditions these models can be reduced to the 
cable model. Continuing with the asymptotic calculations above, we show that 
the cable model can be obtained as an asymptotic limit under assumptions. We 
shall see that there is a hierarchy of electrophysiology models, the Poisson or 
electroneutral models being the most detailed, and the traditional cable model 
being the simplest. 



2 Poisson Model 

We first present the Poisson model, which is essentially equivalent to the model 
proposed in [13] . We consider biological tissue to be a three-dimensional space 
partitioned into the intracellular and extracellular spaces by the membrane. Let 
the biological tissue of interest be divided into subregions fl^ k \ indexed by k. 
We denote the membrane separating the regions fiW and QW by Y^ kV > (Figure 

ED. 

In 0( fc ), the equations to be satisfied by the ionic concentration Cj and the 
electrostatic potential <j> are the following. 

dci 

= — V • fj (ion conservation) (1) 

fi = -A | Vcj + t^V0 j (drift-diffusion flux) (2) 

A<p = — I po + y qzjCi I (Poisson equation) (3) 



-- I po + > „ qzjCj j 
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Here, denotes the flux of the i-th ion. is expressed as a sum of two terms, the 
diffusion term and the drift term. Di is the diffusion coefficient of the i-th ion, 
qzi is the amount of charge on the i-th ion, where q is the elementary charge, i.e., 
the charge on a proton. qziDij (fc^T) is the mobility of the ion species (Einstein 
relation) where Ub is the Boltzmann constant, and T the absolute temperature. 
Fixed background charge density (if any) is given by po, and e is the dielectric 
constant of the electrolyte solution. We note that the above system of equations 
has been used extensively in semiconductor device modeling |22[ [S] and ionic 
channel modeling fiE [ \U \ \W\ . 

Biological membranes consists largely of a lipid bilayer that acts as a capac- 
itor impermeable to ions. In this lipid bilayer are embedded ionic channels and 
transporters through which certain ionic species may pass. With this picture 
in mind, we write down the boundary conditions for the above system to be 
satisfied at both faces of the membrane. 

Consider the boundary condition for the Poisson equation. The value of 
the electrostatic potential and the normal component of the electric displace- 
ment vector D = eE, where e is the dielectric constant and E is the electric 
field, should be continuous at the interface between the cell membrane and the 
electrolyte solution. Therefore, at this interface, 

0(mcm) = ^(Jfc) ( 4 ) 
^(mcm) Q^k) 



(5) 



where 0( mem ) is the electrostatic potential within the membrane, e m the dielec- 
tric constant of the cell membrane, and nA fe " the unit normal at the membrane- 
electrolyte interface pointing from QS k > into the membrane. 

We note that j5]) is not satisfied at the mouths of ion channels. If ion 
channels mouths do not occupy a significant amount of membrane area, the 
above boundary condition may be deemed reasonable. Fortunately, ion channels 
are sparsely distributed even at their peak documented densities }12) . 

The membrane thickness <i m (~ lOnm) is small compared to the curvature 
radius of the membrane and the typical length scale of the system. This implies 
that (j) mom varies linearly as one traverses the membrane from fiw to fi"'. Thus, 

We obtain the following boundary condition, 

where cM fc '- ) = <p^ k ' — <f>"> , C*^ = ^ and n^ fc ^ is the unit normal on the membrane 
pointing from Q,( k > to fl^ 1 '. may be considered the intrinsic capacitance of the 
membrane, which is to be distinguished from the effective membrane capacitance 
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C m to appear later. The jump in the electrostatic potential (j)^ is termed the 
membrane potential and is one of the primary biophysical quantities of interest. 
The boundary conditions for the drift diffusion equations are simple: 

^•11^)= f (8) 

where are ion channel currents carried by the i-th species of ion. We 
note that j\ kl ^ = —jf . These currents can in general be functions of the 
ionic concentrations of arbitrary species on either side of the membrane, the 
membrane potential <// w ) and gating variables which describe the internal states 
of a given ionic channel [TTJ [15] . 

We shall refer to equations H])-© supplemented with boundary conditions 
(J7J) and ©, as the Poisson model. 

3 Non-Dimensionalization and Multiple Spatiotem- 
poral Scales 

We non-dimensionalize the Poisson model. We first rescale the ionic concentra- 
tions Ci and the electrostatic potential <p as follows. 

(9) 

ji = -yqcoji (10) 

where Co ~ 100mmol/l is the characteristic concentration and jqco is the charac- 
teristic magnitude of the transmembrane current per unit area. ksT/q w 25mV 
is the natural unit for the membrane potential. The constant 7 has units of 
velocity =lcngth/timc and its typical physiological range is: 

7 « l(r 5 ~ 10~Vm/ m sec (11) 

We determine a typical length scale of the system. We take equation ([5]) and 
integrate over the membrane surface dQ( k >. 

[ jjidA= [ z i i i -xS k ^dA= f ZiV-iidV 

= - I ZiV-DiCi(ViH)dV 

where we have used dimcnsionlcss variables for ionic concentration and the 
electrostatic potential. In the above, dV and dA denote volume and surface 
integrals respectively and /ii is the chemical potential In Ci + Zi<&. We have used 
the divergence theorem in the second equality and the flux expression ([2]) in 
the third. Let Lq be the typical length over which the flux and the chemical 



<f> = <P, ^ = c Ci 

Q 

Po = QCoPO, £i = Cofj, 



(12) 
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potential vary. Balancing the order of magnitude of the surface and volume 
integrals above, 



Li, 



J {) 



where |9f2( fe )| is the surface area of the region and \fl^\ is the volume of 
We therefore set: 



"-ft-*- -iSr 

The constant Dq k \[imjmsec 2 is the typical diffusion coefficient for ions. The 
quantity I is a measure of the volume per unit surface area, and is a represen- 
tative length scale of the distance between membranes. For a cylindrical axon, 
I corresponds roughly to the diameter of the axon. As we shall see in Section 
18.31 Lo is what is termed the electrotonic length in cable theory. Notice that 
Lq is proportional to This is in agreement with the observation in cable 
theory that the electrotonic length scales with the square root the diameter of 
a cylindrical cable [TT] . 

Given To, we can define a typical time scale To as To = Lq/Dq = l/j. 
This expression tells us that To is cquivalently the time scale in which the 
dimensionless ionic concentration experiences changes of 0(1). We shall call Tq 
the diffusion time scale or the slow diffusion time scale. 

Using To and T) as the representative spatiotemporal scales, we introduce 
the following dimensionless variables. 

x = L X, t = T Q T D , D l = D D l (15) 

i = ^ a = ± (16) 

LiQ J->0 

We can now write the Poisson model CE}-© and (0-© in dimensionless form: 

P = -V x -F, (17) 

OT D 

F, - -A(VxC< + z*QV x $) (18) 

N 

/3 2 A X $ = -(A) + X>Ci) (19) 

i=l 

The boundary conditions are, 

r* w =f , jzm (20) 

z,F, • n^' = al (21) 

Note that a is the dimensionless magnitude of the transmembrane currents 
as well as the dimensionless volume to surface ratio. We have introduced the 
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dimensionless parameters (3 and 9*. The parameter (3 is the ratio between the 
Debye length r d [23| and Lq\ 



P=T, r d=\— — (22) 
L V q 2 c 



The Debye length is typically r d ~ lnm, and is considerably smaller than the 
typical length scale L . The parameter 9* is defined as follows: 

9* = ^ = 9^MEll ~ 10-2 (23) 
a I r d qc r d 

We have, thus, three constants f3, a and 9* that characterize the system. 

Given typical values of I and 7, we can find typical physiological values of 
the parameters [3 and a (the magnitude of 9* is given in (|23| .). Recall that 
I is the (dimensional) volume to surface ratio, and thus, roughly measures the 
separation distance of membranes. Values typical in the central nervous system 
can range from lOOnm to 10/xm. Combining this with the radius of 7 we 
obtain the following physiological ranges for the above parameters. 



L = J— = 10 Mm ~ 1mm (24) 



P = r d\rk- = 10-6 ~ 10 " 4 > Q ' = \/tt = 10-3 ~ 1(rl ( 25 ) 

V lL, V L>v 

We note that while the magnitude of (3 and a depend on the geometry (I) and 
electrophysiological properties (7) of the physiological system under consider- 
ation, 9* defined in (f2"3"|) is a constant that varies little between physiological 
systems. 

We shall exploit the smallncss of the parameter (3 to reduce the Poisson 
model. Note that f3 2 multiplies the Laplacian in (fl8|) . By formally taking 
(3 — > in (|18p . we see that the electroneutrality condition: 

N 

p + z i°i = (26) 

i=l 

should be approximately satisfied in the bulk of the region of interest. The 
electroneutrality condition above is in general not compatible with the mixed 
(Robin) boundary condition of (|20|) . and thus, we have a singular perturbation 
problem which gives rise to a boundary layer at the membrane. Given that 
[3 2 multiplies a second spatial derivative in (JTHJ) , a layer of 0{(3) develops at 
the membrane, where electric charge may accumulate. In dimensional terms, 
this layer has thickness r d ~ lnm near the membrane. We shall refer to this 
layer as the space charge layer or Debye layer. This is a layer that we have no 
need to resolve as long as we are interested in electrophysiology at the cellular 
or subcellular level and not at the molecular level. The biophysical equivalent 
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of this layer in the cable model is the charge associated with the membrane 
modeled as a capacitor, and accordingly, has no spatial extent. Starting with 
Section we shall perform a matched asymptotic calculation that addresses 
the presence of this layer. 

We can now interpret the dimensionless parameter 9* in (|2U|) as follows. The 
constant ksT/q is the typical magnitude of the membrane potential, whereas 
qcQTd is a natural unit of surface charge density since gives the surface charge 
thickness. Thus, ^^/(fc^T/g) is a natural unit of capacitance per unit area. 
The constant 9* expresses the membrane capacitance per unit area in these 
natural units. 

Before we can perform asymptotics on the model, we would like to identify 
other spatiotemporal scales that the Poisson model possesses. Differentiate both 
sides of equation (fl9|) in td and take the integral over fiW. The left hand side 
yields: 

m ot d J aQ(k) 6t d V dn hl J 

(39*— dA. 

We used the boundary condition (|20p in the second inequality. The right hand 
side yields: 




r N 
dV= aS^3i dA. (28) 

Janw ~[ 



where we have used p7|) . (|20p and the divergence theorem. The above says that 
the change in total charge within QP*> comes from transmembrane currents. 
Balancing the quantities in (|2"T|) and (f2"B")) . we see that the membrane potential 
and hence the electrostatic potential can vary on the time scale of f3^T$. It 
is an interesting coincidence that 9* and a are roughly of the same order of 
magnitude, as can be seen from ([23]) and (|25| . Thus, this time scale is roughly 
equal to (3T$, which we shall call the membrane potential time scale. Given 
the smallness of /3, the membrane potential time scale is considerably smaller 
than the slow diffusion time scale To. We shall see in Section 18.31 that the 
membrane potential time scale {3Tq corresponds to the "diffusion" time scale of 
the membrane potential in the traditional cable model. 

There is yet another time scale, which corresponds to charge relaxation: 

d ( N \ N 



ZiCij + other terms 

(29) 



where we have used the Poisson equation (fT§|) in the last equality to replace 
A$. We see that charge density decays exponentially with a time constant of 
(3 2 Tq = t^/Dq = lnsec. We can infer that this time scale is only important 
where the electrolyte solution may deviate significantly from electroncutrality, 
i.e., within the space charge layer. 

We thus see that there are three time scales present in the Poisson model, 
To, /3Tq and P 2 Tq. We list the physiological values for these time scales. 

T = 1CT 1 - 10 3 sec, (3T = 1CT 2 - 1 msec, /3 2 T = 1 nscc (30) 

The time scale of greatest interest is the (3Tq time scale, in which the membrane 
potential varies. This is also roughly equal to the time scale in which the 
most rapid physiological processes take place, such as channel gating, chemical 
neurotransmission and calcium concentration changes [I]. We shall thus focus 
our attention on this time scale and rescale the time variable td to a newly 
rescaled time variable ry = td/(3. We write Ci,3> as functions of ry rather 
than td- Equation (TIT]) is rescaled to: 



The (3 2 To time scale and the space charge layer within which this time scale is 
relevant are spatiotemporal details that we have no need to resolve. The To time 
scale is important with regard to long term changes in ionic concentrations. We 
shall make some brief remarks about this time scale in the final section. 

An overarching goal is to computationally investigate the three dimensional 
electrical activity of complex physiological systems. A great difficulty with 
the Poisson model is that one inevitably needs to resolve spatiotemporal scales 
associated with the space charge layer in a numerical simulation, making such 
computations prohibitively expensive. It would therefore be computationally 
desirable to have a model that resolves the membrane potential time scale but 
does not resolve the Debye spatiotemporal scales. 

4 Electroneutral Model 

We propose the following as a computationally efficient alternative to the Pois- 
son Model: 



dry 



/3V X • Fi 



(31) 



= ^i+/3Vx-F i 

OTy 

\ = -ACVxCi + ZjCiVx*) 



(32) 



(33) 



N 




(34) 



z i F l ■ n 



(ki) 




+ Oiji 



(35) 



dry 
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The Poisson equation in the Poisson model has been replaced by the electroneu- 
trality condition (|34[) . Since this is an algebraic condition, it does not require a 
boundary condition at the membrane. The boundary conditions for the drift- 
diffusion equations (|32|) and (|3"3"|) are given by ([33)) . In comparison to Pl"j) . we 
have an additional term: 

(Ti is the amount of electric charge at the membrane face contributed by the i-th 
species of ion. In the electroneutral model, the electric charge within the Debyc 
layer is represented as a surface charge density of zero thickness. In this picture, 
the amount of ionic current gzsfj • n either contributes to the change in surface 
charge density Oi or flows across the membrane through ion channels. This 
picture is better aligned with the biophysical view of the membrane in the cable 
model, in which the membrane is a capacitor within an ohmic medium. One 
important advantage of the boundary condition ([35]) compared with ([2Tj) is that 
the parameter values in (|35[) are directly observable experimentally. Since the 
Debye layers are too thin to be explored experimentally, the parameter values 
in (p?Tj) can only be inferred, as argued in detail in [TS] . 

The surface charge contributions a% must be related to the dynamic variables 
Ci and/or to close the system of equations. First we let 



N 



7 = 9<5> (kl> . (37) 

i=l 



This relation says that the total amount of surface charge a is linearly propor- 
tional to the membrane potential where 9 is the effective dimensionless 
membrane capacitance. Note that 9 is different from 9* , the intrinsic membrane 
capacitance, used in ([20"]) . The lipid bilayer sandwiched by the two boundary 
layers considered as a whole gives rise to a capacitor with the effective capac- 
itance 9. This is the capacitance that is measured experimentally, given that 
it is impossible to to distinguish the contributions to the capacitance from the 
Debye layers and the lipid bilayer. The relation between these two quantities 
will be clarified in Appendix [TT] Now, define A 4 as the fraction of the total 
charge a that is contributed by the i-th species of ion: 

CTj = \icr. (38) 

We let Xi evolve according to the following: 

9X t _ Xi — A 4 ~ zfCi , 

^ = — Xl = ^L4^ (39) 

The charge fraction Aj relaxes to A^ in the charge relaxation time scale. The 
specific form of A 8 was derived in [15j , but is also given in Appendix 1111 Note 
that: 

N \ N / N \ 



d 




^(IJM =D A *- A *) = 1 - 
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and therefore, $2 i=1 K = 1 provided that 5Z i=1 Aj = 1 at the initial time, as 
required by the definition of Xi as the charge fraction. In [15j . Xi was used in 
place of Xi in ([55)) , in which case the charge fraction relaxation equation in (|39|) 
is not needed. This original system, however, leads to ill-posed behavior which 
we examine in Appendix 1121 

We shall call the system (f3"2" j) - (|35p and (|39[) the electro-neutral model There 
is no longer a space charge layer to be resolved at the membrane, since the 
presence of the surface charge has been taken care of in the boundary condition 
(|35|) . The charge relaxation time scale only appears in a simple ODE (f39|) . and 
does not pose serious difficulties in the construction of a numerical scheme [16] . 
We propose the electroneutral model as a computationally tractable model that 
addresses the shortcomings of the cable model pointed out in Section [TJ 

An important difference between the electroneutral model and the Poisson 
model is what the state variables are. In the Poisson model, specifying the ionic 
concentrations at every point in space is enough to describe the state of the 
system. The electrostatic potential can be found from the ionic concentration 
profile by solving the Poisson equation (fl9|) with the boundary conditions (f20|) . 
The difficulty, though, is that we must specify the ionic concentrations up to 
the boundary to within the space charge layer. The electroneutral model, on the 
other hand, does not require the ionic concentration profiles in the space charge 
layer. The spatiotemporal details of the space charge layer are represented 
by the the membrane potential $( fc ') and the charge fractions A^. The state 
variables for the electroneutral model thus include the ionic concentration profile 
as well as the membrane potential $( fci ) and the membrane charge fractions Ai. 
This means in particular that we need to specify the values of these quantities 
as initial conditions. 

In the electroneutral model we have ion conservation in the following sense: 

^-([ z l C l dV+ [ peX^&^dA) =- [ pajidA. (41) 

This equation says that for each ionic species the change in the sum of the 
ionic content of the region £l' fe ) and of the space charge layer is equal to the 
transmembrane current that flows out of this region. This is an important 
property not only from a physical point of view, but also from a practical point of 
view if we are to perform long-time calculations of ionic concentration dynamics. 

The natural question that arises is whether the electroneutral model is in 
any way an approximation to the Poisson model. We investigate this question 
using both asymptotic and numerical computations. Beginning with the next 
section, we present a matched asymptotic study to show that the electroneutral 
model gives an approximation to the Poisson model. In Section [9l we shall 
computationally investigate how well the electroneutral model approximates 
the Poisson model. 
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5 Matched Asymptotics 



We recall the Poisson Model: 



dr v 
F 



A(V x Ci+^aV x $) 



/3V X • F, 



(42) 



(43) 



Po + X] ZiCi 



(44) 



i=l 



Recall from (|3"Tj) that we rescaled time to Ty to capture the dynamics in the 
membrane potential time scale. We now perform matched asymptotics on the 
above to clarify the relation between the electroneutral and Poisson models. 

As noted earlier, a boundary layer of thickness 0{f3) develops at the mem- 
brane when /3<1. We therefore introduce an inner layer of thickness 0{(3) at 
the membrane. We shall continue to use the terms space charge layer or Debyc 
layer to denote this layer. 

We need to introduce another spatial scale of order 0(^/]3) at the membrane. 
This need arises as the result of introducing a newly rescaled time variable Ty. 
The spatial scale of order 0(y/]3) corresponds to the length over which ions 
can diffuse in the membrane potential time scale, [3Tq. Formally, the necessity 
for this layer can be seen by noting that (3 multiplies a second order spatial 
derivative in (|3~Tj) since F; is itself written in terms of spatial derivatives (c.f. [T8"|) . 
We shall refer to this layer as the intermediate layer or the fast diffusion layer. 
It is interesting to note that the presence of such layers have been postulated to 
account for K + ion accumulation in the extracellular space of the squid giant 
axon [3] . We thus have three regions to consider in the asymptotic calculations 
to follow: the inner and intermediate layers located adjacent to the membrane, 
and the region away from the membrane, which we shall call the outer layer. 
We perform two matching procedures, at the inner-intermediate layer interface 
and at the intermediate-outer layer interface. We have summarized the relevant 
spatial scales in Figure [2] 

The above discussion prompts us to expand the physical variables in powers 
of \f]3 instead of (3: 



The other two parameters of the system, a and 9* are also small (c.f. (|25l) . (|2"3|0 . 
but we shall treat them as being 0(1) with respect to (3. We note that (3 is 
typically a few orders of magnitude smaller than a or 9*. The smallness of a 
and 9* will be later exploited, in sections 18.31 and [6] respectively. 

In performing matched asymptotics at the membrane, we introduce a coor- 
dinate system at the membrane £ = (£i, £2, £3), where the £1 axis is taken to 
be perpendicular to the membrane, while £2 and £3 are curvilinear coordinates 



Ci(X,ry) 

*(X,7V) 



C°(X, Ty) + y/0C}(X, Ty) + /?C?(X, Ty)--- 
$°(X, Ty) + vW(X, Ty) + /3$ 2 (X, Ty)--- 



(45) 

(46) 
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Figure 2: A schematic of the relevant spatial scales used in the asymptotic 
calculations. The solid lines denote the membrane and the dotted lines are the 
interfaces between the boundary layers. The inner-most layer has width f3, the 
intermediate layer y/]3. The typical membrane separation is a and the typical 
length scale associated with the membrane is l m . £ is the membrane- fitted 
coordinate used in the matched asymptotics calculations. 

that run "parallel" to the membrane. The £1 axis will be rescaled to yield coor- 
dinates in the intermediate layer £ a such that £1 = v^Ci and in the inner layer 
£ 6 such that £1 = 

We must now ask how we are to rescale £2 and £3. There are at least two 
spatial scales that are relevant: p K the dimensionless curvature radius of the 
membrane and lj the dimensionless length scale on which one may see 0(1) 
changes in ion channel current density. Let l m be the smaller of the two spatial 
scales lj and p K . We shall call l m the membrane length scale. The question 
raised at the beginning of this paragraph can be answered by comparing the 
relative magnitude of this length scale to the 0(\fj3) length scale. 

If l m is considerably larger than \fj3, there is no need to rescale £2 and 
£3. If l m is order 0{\ff3), we must scale £2, £3 to £2 = v7^2' & '£3 = V/^' 6 so 
that the curvature correction and the ionic fluxes parallel to the membrane are 
O(l) quantities when measured in the intermediate layer coordinate £ a . We 
shall mainly be concerned with the case l m > \/~j3 but we shall quote results of 
calculations when l m ~ yfft. 

We point out that there could be situations in which l m is small only along 
a certain coordinate direction. For example, if we take a cylindrical axon with 
diameter 0{\fj3), and take £2 to be the angular coordinate, and £3 to be the 
axial coordinate, the curvature radius along the £2 coordinate is 0(^/j3) whereas 
the curvature radius along the £3 coordinate is large (curvature is negligible). In 
such cases (and if the cylindrical axon is endowed with near uniform ion channel 
density so that lj is large), we need only rescale £2 but not £3. We shall not 
deal with such cases, since such an analysis will follow along similar lines as the 
case in which l m ~ yffi. 
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6 Inner-Intermediate Matching 



Wc first consider inner-intermediate matching when l m > yfft. 

Consider the membrane surface facing fl^ k ' . We now introduce a coordinate 
system £ so that the £1 coordinate direction is perpendicular to the membrane. 
We let £i=0 coincide with the membrane face, and let the positive £1 axis 
point into the region For simplicity, we shall assume that the membrane is 
flat, i.e., that it has no curvature. Therefore, we can take the coordinate system 
to £ to be orthonormal. When l m > i/]3, it turns out that curvature corrections 
produce only higher order terms that we can ignore. 

In the inner layer, we rescalc £ as: 



F b 



£3 



(47) 



The equations are: 



dry 



-A: 



oc b 



dF b 



P 2 



Po 



N 

E 



Co 
. i 



2 



9F% 

oe 3 



(48) 

1,2,3 (49) 
(50) 



Since the inner layer is adjacent to the membrane, we must supplement the 
above with boundary conditions, suitably rescaled: 



$(0) 



z-F b I 



(3a ji. 



(51) 
(52) 



We shall make the simplifying assumption that the transmembrane ionic current 
densities ji are given functions of position (on the membrane) and time instead 
of being functions of Cj, $( fc ^ and the gating variables. 
In the intermediate layer we rescale £ as: 



The equations are: 



tx = y/P$, 6=8, 



dry 



Fi, 



-D 



dFl 



dF? 2 



dtf 9£ 3 a 



(53) 



(54) 



<9 2 <f> a 



■0 



2 / d 2 $ a 



p 

A' 



i=l 



p= 1,2,3 (55) 
(56) 
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Substitute (|43|) and (|4*S|) in the inner layer equations (|4*8"|) - (f5"01) . and collect 
like terms in order (3. The expansions of Cj and 4> in \f]3 induce expansions of 
¥i in terms of \f]3. We shall denote the 0{y/~]3 ) term as F*\ For example, 

^ = _ A ^ + „ Cf ^ + , iCr ^ (58, 



By applying the same procedure to the equations (|54[) - (|56p . we may obtain 
analogous expressions in the intermediate layer. 

We derive matching conditions at the inner-intermediate layer interface in 
terms of the ionic fluxes. Note from (|48|) and ([52]) that: 

(59) 
(60) 



From this, we find that 





= 0, 


fa 


fj=0 


= 


dF$ 






e?=o 


= 




7716O 

Fl — 




E 





(61) 



within the inner layer. 

Now, consider the p = 1 component of (|4"9"|) and (f55j) . i 7 ^ and F^ . We intro- 
duce a matching coordinate system (J* in between the inner and intermediate 
layers such that, 

= n(B)C, lim ^ = lim 77 = (62) 
/3— >o rj p^o 

Applying Kaplun's matching condition [SUTU] to and F^, we obtain: 

W» = (63) 
ft (^ F ™ (> ') + * (^ $ ") " F " (,f °) = ,64) 

--L^o^-FS 1 ^")) =0 (65) 

Condition (|63|) is automatically satisfied by ([61) . Condition (|64|) . taken together 
with (fBTj) . yields: 

hmi^H = ^f L- = 0- (66) 

/3^0 'Si— u 

We thus have the matching condition for the leading order ionic flux in the 
intermediate layer. The last matching condition (|65[) . combined with (|6ip . yields 
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the following. 



Urn (F<f \^L? >) ^FtiivC) - F.a\vC)) = 0. (67) 

To evaluate (f6"T|) , we need to calculate Ci and $ to leading order in the inner 
layer. From (fSTj) . ([50]) and (|5ip we see that the leading order terms satisfy the 
following one dimensional boundary value problem in in the inner layer: 

60 



°~^T + zA ^T (68) 



2^60 



a 2 $ 



AT 

(pb + $>Cf) (69) 

i=l 



* M (d: = o) - $ (,),60 (^ ,(0 = o) 



as 



(70) 



C*°(tf = oo)=Cf(£=0) (71) 

$ bo (ej = oo) = $ a0 (ei =o). (72) 

The last two conditions come from matching conditions at the inner-inter- 
mediate layer interface. Here, £i refers to the inner layer coordinate system 
on the side of the membrane r^ fc ^(note that we are now working on the f^* 1 
side). Equations (f6"5)) . with boundary conditions ([70)) - (|72p can be solved 
explicitly under the approximation that 9* is small, a reasonable approxima- 
tion since 9* rts 10~ 2 (cf. (|2"5|) V We quote the results below, and relegate the 
calculations to Appendix [TT1 

= * o0 (0)-£«p(-rtf) (73) 

Cf te b ) = Cf (0) (l + ^ exp(-rtf)) (74) 

r 2 = ^z 2 cf(0), r>o (75) 

5* = *'^ 2 w <x = \a (76) 
- e$ (fci),ao I _ J_ i _L + J_ (77) 

9* 1 j 

In the above, a denotes the total charge in the Debye layer, and a% is the 
charge contributed by the z-th species of ion. Thus, &i/<r = A; is the charge 
fraction contributed by the z-th species of ion. Note by design that Yli=i = 
1. The variable 9 is the dimensionless effective membrane capacitance to be 
distinguished from the dimensionless intrinsic membrane capacitance 9*. Wc 
refer the reader to Appendix [TT] for further elaboration. 
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Now, consider Q48p and ([5^)1 at the first non-trivial order: 



«i^?U=o = ai< (78) 



dry ' d& ' 1 11 

Since our goal is to evaluate (|67|). we would like to obtain an expression for F^. 
We integrate the above in to obtain: 

-z,Fg = aj, + z, f' 

. ali + ,220^ + f' JL^mi^-T^ (79) 
ot v J 3t v r 

^ a j j;+2i ^!M^ + /charge 
OTV 

where we used (fT4"f for C*°. We can finally consider condition (pT)) . We would 
like (jUT)) be satisfied regardless of the value of £ v . For the i 7 ^ 2 term, taking 
j3 — > in ([57| amounts to studying the behavior of (f?T)|) in the limit £j ~~ * 00 • 
Take $ ->■ oo in / chai -g c . 

lim Zi I chaise = (80) 

£*>->oo C/TV 

where we have used (|75|) . We thus conclude using ([75)1 and the above that: 



ZiF% = ( P- + aj« ) + Zi^p^J + 0(exp(-r^)) (81) 



a* , . V 9Cf°(o) eb 

5 1- aji + Zi 

OTv J OT\ 

Note that Fff thus consists of a constant and a linear component in £j as well 
as a residual term that decays exponentially. We now expand the intermediate 
layer expressions of (|57|) at = 0. 

_ nc ^f(o) , „ x 



(0) 



where we have used to eliminate F^f(0). Substituting the above as well as 
(Ell into (1571). 



(83) 



lim fe(0) + i (|^ + aj, 

/3^0 V Zi \OTv 
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The necessary conditions for the above to be satisfied are 



^(0) 
dCf(0) 




d&i 



(84) 



dry 



d*g°(0) 



(85) 



The second expression ([55]) is automatically satisfied as can be seen by tak- 
ing ([55]) to leading order. Equation ([84]) together with ([66]) are the matching 
condition we set out to obtain. 

When l m ~ v^9, as discussed at the end of the previous section, we must 
rescale the coordinates so that £ p = \//3£p ', p = 2,3. We can obtain the 
matching conditions for this case in a manner similar to the l m > ^/J3 case, 
although the calculations are more involved. The matching conditions corre- 
sponding to ([84]) and ([85]) are respectively [16] : 



Here, the operators Vs» and Vs<»- denote respectively the gradient and diver- 
gence operators on the membrane, where the length is measured in terms of 



Compared with ([84]) . equation (|86|) has an additional membrane drift dif- 
fusion term. The surface gradient of the chemical potential potential ^ = 
InCf + Zi$ a0 , scaled by the diffusion coefficient, gives the drift velocity of ai 
along the membrane. 

We shall henceforth limit our attention to the case l m > (3. 

7 Electroneutral Model as Approximation to Pois- 
son Model 

We now examine the relationship between the Poisson and electroneutral mod- 
els. Consider two pairs of ionic concentrations and electrostatic potential Cf , <& EN 
and Cf°, $ Po , which evolve according to the electroneutral and Poisson models 
respectively. Wc postulate an expansion of Cf N ,Cf° and $ EN ,<i> Po in a/9 of 
the form ([45]) and (|46[) respectively and see if the electroneutral and Poisson 
models produce the same leading order equations. 

First consider the intermediate layer. We write equations ([32l)- (|35[) and (1391) 
of the electroneutral model in the £ a coordinate and write out the leading order 




(86) 



(87) 
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equations. The 0(1) equations are: 



dry 8Q 



(88) 



p aO,EN _ f, I I ^aO,EN 



1 dH J 

N 

= /J + ]>>C? ' EN (90) 
i=i 

^° lEN (^ = 0) = (91) 



il 

The O(v^) equations are: 

8Cf ' EN _ 9^ a 1 1 ' BN 



dry flfj 



(92) 



o^al.EN a fF( aO,EN n^a^EN 



Fal,EN t-v | 
n = -Di I — — h ZiCi — — h 2iC 



(93) 

AT 

= /> o + I>Cf- EN (94) 
jtf^K? = 0) = + otf, (95) 

a\ \EN \ 2/-<aO,EN 

OA i _ A i ~ A i \ EN _ z i °i / Qfi \ 



dry " /3 ' 1 , z 2 ,0f < EN 



The same procedure on the Poisson model yields the following. The 0(1) equa- 
tions are: 



dcf Po _ dF^ Po 



(97) 



dry 

( Q^aO.Po a ( f ( aO,Po , 

^°' P ° = -A + ^^°' P °^ !r j 08) 

N 

= pb + 5>Cf lPo (99) 



J^°.P°(g = 0) = ( 10 Q) 

Equation (|100p comes from the matching condition (|66|) . The (D(\fj3) equations 
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are: 



dCf Po _ dF^° 

, ' r)r al ' P ° «<T,aO,Po 




N 

Po 



o = p a + J2^cf 



i=l 



Ki ' (e? = 0) = + aji 



A7° = 



z. 



2^iaO,Po 




where equations (|104|1 and (|105[) come from the matching condition (|84[) . 

We sec that dHHJ) - dHU) , dM))-® arc identical to ([W l - flUO)) . (fTUTjl - lfTTHj) . ex- 
cept for the difference between A^ and Aj in equation (P5|) and (|104|) . In Ap- 
pendix we show that in fact (Eq. (|205[) ): 

^L_ = ^i_ +0(/3 ) (106 ) 
ary dry 

Therefore, Af N may be replaced by Af N without affecting the order of the ap- 
proximation. This shows that Cf™, < & EN and Cf°, $ Po satisfy identical equations 
in the intermediate layer to order 0(^/]3). 

The same procedure in the outer layer shows that the two models agree 
up to equations of order 0(/3 3//2 ). We thus see that the electroneutral model 
formally approximates the Poisson model in the intermediate layer and outer 
layers, where the biophysical processes of interest take place. In Section [9l we 
shall show computationally that the electroneutral model indeed provides an 
excellent approximation to the Poisson model. 



8 Equations in the Outer Layer 

We continue with the asymptotic calculations with the goal of obtaining the 
cable model under certain conditions to be set forth below. 

8.1 3D Cable Model 

We now consider intermediate-outer matching. We first turn to the equations 
satisfied in the outer layer, which can be obtained by substituting (|4"1)|) and (|4"5j) 
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into (5U), CEED and (32). 



S^.o, H_o (107) 

-V X -F° (108) 



c?7y 

iV JV 

o = pb + Yl ziC i = Z z * c * 2 ( 109 ) 



i=l 



Equation (|107|l tells us that Cj to leading order does not change in the ry time 
variable. We still need to know the evolution of $°. This can be obtained by 
summing (| 1 08(1 in i and and using (|109|l to conclude: 



(110) 

=1 / 

This is the equation satisfied by <3?° in the outer layer. In order to obtain the 
boundary condition for this equation, all we need is (Y^iLi z i^i) ' n^ fc ^. 
Let J = J2"=i z i F i- We sha11 

use the usual subscripts and superscripts on 
J to denote terms of the expansion of J in (3 in the different layers, induced 
by the expansion of Fj. We find from (|54")l and ([51?)) that: 

N N 

5>Cf = 0, 5><7f = (111) 

(112) 



i=i »=i 

dCf _ dFg dCf dF^ 1 



dr v " dQ ' dr v dQ 

Using the above relations we see that: 

djf _ dj{ 



= (113) 



The value of Jf° and J£ x at = can be computed from (|55J) and l[8"4"|l. and 
we see from (| 1 13[) that: 

as,{H),aO N 

Jf = 0, -Jf = 6 — + £ aj< (114) 

where we have used y ■_ 1 A; = 1. Following the same matching procedure as 
for the inner-intermediate layer matching, we conclude: 

Q$(kl),aB 

^ + «E^ = ^ ( 115 ) 
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We can now use the above as the boundary condition for (| 1 1 0[) and explicitly 
write down the equations satisfied in the outer layer. 



V • (AV$° + V-B) = 



(116) 



(AV$° + VB) ■ n< w > = 6— + ah 



(117) 



OTv 

N N 



N 



ion — 



^2 ji 



(118) 



i=l i=l 



i=l 



Note here that A and B are functions of X only, and do not depend on time, 
since C° does not change in the Ty time scale. 

There is one difficulty here that needs to be pointed out. Equation (| 1 1 5(1 
and (|117[) are not exactly the same. In (|115p . is evaluated just outside the 
inner layer, whereas in (jTTTJ) , $( fc is evaluated just outside the intermediate 
layer. There is a similar concern for the transmembrane current terms ji if they 
are functions of Ci or $(*0 . 

From (jllip and (|6"6")) . and the fact that Cf° and 3> a0 must match to leading 
order at Q = oo to the outer layer solution, we see that Cf° and <E> a0 decay to a 
uniform state after an initial transient (note Cf° decays to a constant where as 
i> a0 decays to a time- varying uniform state, whose value is equal to <I> (£j = 0)). 
Therefore, after an initial transient, the discrepancy between <J>( H )'°, $( fe 0. at) an( j 
Cf, C° will decay to 0. 

This model is valid to leading order outside the intermediate layer of thick- 
ness 0(y/~j3). We shall call this the 3D-cable model 

The the 3D-cablc model may be derived very easily from the electroneutral 
model. Consider equations (|32[) - ([551) of the electroneutral model. We can take 
the time derivative of the electroncutrality condition and substitute (j3"2")) to 
obtain the elliptic equation satisfied by the electrostatic potential, (|116p . Sum 
(|35|) in i and we find the boundary condition (|117p . The coefficients A and B 
in equation (|116p are now time dependent, but we can see from (|32p . that to 
leading order, the ionic concentrations do not change in the membrane potential 
time scale. Thus, A and B are constant to leading order. The ease with which 
one can see the correspondence between the electroneutral model and the cable 
model is an appealing feature of the electroneutral approach. 



We reach a further simplification by considering the following situation. Suppose 
the long time average of the transmembrane currents j% is equal to 0. That is 
to say, if we average over a sufficient long time, there is no net current flowing 
through the membrane. An electrically active cell whose ion channel currents are 
quickly counter-balanced by ionic pumps may fit this category. Then, the ionic 
concentrations should relax to a stationary value in the slow diffusion time scale. 
If there arc no fixed charges p, or if the fixed charges are spatially uniform, the 
resulting ion concentration profile will be spatially uniform within each region. 



8.2 Simplified 3D-Cable Model 
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We apply the above 3-D cable model to this situation. From (|116[) - (|118[) . we 
obtain: 

A$° = (119) 

f)(f)(kl),0 

-A™ V$° ■ n< w > = 6— + aI ion (120) 

OTy 

The gradient of B vanishes because of the spatial uniformity of Cf. Note that 
is a constant that depends only on the region number (k), and expresses the 
ohmic conductivity of the electrolyte medium. We shall call this the simplified 
3-D cable model. We note that this system, when homogenized in a quasi- 
periodic domain, gives rise to the bidomain equations, which are widely used in 
simulations of organ- level cardiac electrophysiology [17l QT] . 



8.3 Derivation of Standard Cable Model 

We now derive the traditional cable model by considering the above simplified 
3D-cable model under specialized geometry. We note of an analysis of a similar 
situation for a passive cable in which a different approach is used to address this 
issue [2"T] . 

Consider an infinitely long cylinder of radius r mt . The dimcnsionless radius 
will therefore be rj = r mt /Lo- This infinite cylinder is surrounded by an ex- 
tracellular space which lies between this cylinder and a concentric cylinder of 
radius r cxt (> r mt ). This extracellular region is insulated at the outer boundary. 
We shall let £ = r oxt /r mt . Equations of the simplified 3D-cable model (|119[) and 
(1120)) specialized to this situation are: 

in fi int ,^ GXt (121) 
0§£l + a I ion , [$] = $ int - $ cxt at R = r/* (122) 

OTy 

at R = 1^ (123) 

To avoid cluttered notation, we have eliminated the superscript 0. In the above, 
R is the radial, Z the axial coordinate and Ad denotes the Laplacian on the 
plane Z = const. Equation (|122p is satisfied at R = r\ from both the intracellular 
and extracellular (r; + ) sides. The superscript k denotes either the intra or 
extracellular region. 

We shall now take 77 to be the small parameter in our system. What follows 
is a thin-domain asymptotics calculation used for example in lubrication theory 
[7]. We rescale the the radial coordinate to R = r/p in p2ip ~ (|123p so that the 



dZ? 



Ad* = 

A dR- 

A c*td± = 
OR 
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cell membrane corresponds to p = 1. 



i 

Tj dp dry 
d$ 
dp 



V 2 — + A f) $ = Qmn int ,n c * t (124) 
4(*0 r)<T> r>[<61 

oJion, [*] = * - * at p = 1± (125) 

c/y 



where A^, denotes the rescaled Laplacian on Z = const. We now expand <& in 
powers of rf in the following fashion: 

$ = $° + ?; p$ 1 + ... (127) 

We let p = 2 so that we obtain nontrivial expressions when the above substituted 
into (T2H) : 

A^$° = (128) 

a 2 $° 



^tA 5 $i =0 (129) 

Consider the boundary condition (|125[) . Upon substitution of (|127p . we see that 
a distinguished limit can be obtained by taking a ~ 77. This is in fact, hardly 
surprising. In Section [3l we introduced a as the volume to surface ratio of the 
domain of interest. The dimensionlcss radius 77 is exactly equal to this ratio (up 
to a factor of order 1). We shall thus take 77 = a. Therefore, 

a$o,(fc) 

— = 0atp = 1 ± j/9 = £ 130 ) 

dp 

n.d® 1 6 d\$°] d® 1 

-A( fe )^ = -^i+J ion at /9 =l ± , ^L=0at P = e (131) 

dp a dry dp 

First of all, (fT2g|> with ([HQ]) tells us that $° is constant for fixed Z. In order 
to find the Z dependence of $°, we need to look at the next order, Q129p . 
The solvability of this equation with respect to "i" 1 requires that the following 
identities between an area and a line integral hold for each Z = Zq. 

[ A £) $ 1 dA = [ ^-ds (132) 

Jp<l,Z=Z a Jp=l-,Z=Z VP 

f d® 1 f d® 1 

A B <S> 1 dA = - ds+ —ds (133) 

i<p<€,z=z J p= i+,z=z dp Jp=£,z=z dp 

where dA denotes an area integral and ds denotes a line integral. Applying the 
above to p^]) and (fHiTj) we find that: 

^ ^^= 2 ^-^ J +/ J ion# (134) 



dZ 2 a dry 

d 2Q0, C xt d\<f> Q ] 



-7r(^-l)A ext =^-~^+ I 4>n# (135) 
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Dividing by the prefactors and adding the two expressions, we immediately 
obtain the cable equations: 



dZ 2 \a dry 

- m = l- 1 i ion d^ (137) 



G eS nA int 7r(fa _ ' loin ' 2n 

We have thus succeeded in deriving the cable model. We note in particular that 
Z is measured with respect to the length scale Lq, which we can now identify as 
the electrotonic length. The time variable ry is measured with respect to (3Tq 
which tells us that j3T$ is "diffusion" time scale for the membrane potential. 

If 1 <C £ <C a _1 (r lnt <C r cxt <C Lq), we can take the extracellular space 
to be an isopotential compartment and set G eff = 7rA lnt without sacrificing the 
validity of the above cable equations. In dimensional terms, the above equations 
take the following familiar form: 

Pm Cra—KT + ' P m = 2lTr ( 138 ) 



R dz 2 ^ V 9t 

R = R int + i? oxt (139) 

^ = ^1^^ ^ int = ^ int ) 2 (140) 

1— 1 

_1_ = scxt J2 { ^f^cT\ = n((r^ ~ (r int ?) (141) 

2— 1 ^ 

Here, (f> m is the membrane potential and i lon is the dimensional transmembrane 
current, averaged over the Z =const cross-section of the membrane. 

We note that the above derivation of the cable model did not assume an 
axisymmetric solution to the equations. The axisymmetry, or more strongly, the 
constancy of the electrostatic potential for each Z cross-section is a consequence 
purely of the scaling relations. Related to this is the observation that the above 
can be generalized to arbitrary cross-sectional geometry. All we have used is the 
divergence theorem as applied to each cross section; we have made essentially 
no use of the fact that the cross-section was a disc. 



9 Numerical Validation of Asymptotics 

In this section we shall test the behavior of the electroneutral model against that 
of the Poisson model by way of numerical simulations. We confine numerical 
validation to test cases which reduce to one dimensional computations. This 
is because the Poisson model requires extremely small time steps and spatial 
resolution, which makes it computationally overwhelming to compare the two 
models in a full two or three dimensional setting. We have considered two 
geometrical situations, one spherical and one planar, but we shall only discuss 
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the spherical case, since results for the planar calculations are very similar to 
the spherical [16]. 

We take a spherical cell of radius I. Let the center of the cell be the ori- 
gin, and let r be the radial coordinate. We seek solutions to the equations 
(electroncutral or Poisson) which depend only on the radial coordinate r. We 
have thus a one dimensional problem. The region characterized by r < I is the 
intracellular space. We confine our simulation domain to r < 21 and impose no- 
flux boundary conditions at r = 21. Thus, our extracellular space is the region 
/ < r < 21. The I we use here as the radius of the cell is to be identified with 
the I we introduced as the volume to surface ratio in Section [3] 

We now rescale length so that Lo, (HI}, is the representative length scale. 
The dimensionless cell radius is now a = 4-. We shall continue to use r as 
our dimensionless coordinate. Thus, r < a is the intracellular region and a < 
r < 2a is the extracellular region. We use the finite volume method to perform 
the simulations. We subdivide the computational region into spherical shells 
indexed by k. The thickness of the spherical shells is made to be smaller near 
the membranes so as to resolve the space charge layer and the fast diffusion 
layer. The details of the numerical scheme as explained in |16j will be reported 
elsewhere. 

We consider four ionic species with the following dimensionless diffusion 
coefficient and valence. 

Di = 2, Da = 1/2, #3 = 1, £ 4 = 1/2 (142) 

zi = 1, z 2 = 1, z 3 = -1, z 4 = 2 (143) 

Recall from Section [3] that f3, a and 9* are the dimensionless parameters that 
characterize the system of equations. The parameter 6* = 10 -2 (cf (|23p has a 
fixed value. We consider three pairs of parameter values: 

(A a) = (Kr 3 ,ifr 2 ), (ht 3 - 5 , ht 1 - 5 ), (lcrMtr 1 ) (144) 

We expect the electroneutral model to be a good approximation to the Poisson 
model for small values of /?. We thus take (5 to be slightly larger than the 
typical values f3 = 10~ 4 ~ 10~ 6 to perform a more stringent test of validity of 
the electroneutral model. 

We shall start our simulation at time ry = —T r where T r is positive. The 
reason for this will become clear shortly. For the electroneutral model, we set 
the initial conditions at Ty = —T r for Cj to be: 

C 1 (r,-T r ) = l + O g (2\a-r\-a) C 2 (r, -T r ) = 2 - C x (r, -T r ) (145) 
C 3 (r,-T r ) = 2 (146) 
C 4 (r, -T r ) = 10~ 6 for r<a C 4 (r, -T r ) = 1CT 3 for r > a (147) 

4 

Po( k r) = -^2z i C i (x,-T r ) (148) 
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We let C g = 0.9 so that there is a steep initial gradient of the ionic concentra- 
tions. The very small initial values of C4 are motivated by calcium concentration 
profiles in physiological systems. At the membrane boundary, we must specify 
&m{i~v) and At(a±,7v) at t v = —T r . 

$m(-T r ) = $(a-) - *(a+) = ~ (149) 

a 

where a+ and a— denote the r > a and the r < a faces of the membrane 
respectively. 

For the Poisson model, we need to specify the initial ionic concentrations. 
Given initial conditions for the electroneutral model, we set the corresponding 
initial conditions for the Poisson model to be: 

d(r, -T r ) = Cf cct " al (r, -T r ) ~ ^^"T^f ^ ^ r < a (151) 

d(r, -T r ) = V, -T r ) + ^f^ ( t\ 3 Tr)e * iir>a ( 152 ) 

Zi(4:Tr/S)((2ay — a A ) 

The rationale for setting C, as above is the following. The initial conditions 
for the electroneutral model says that each ionic species contributes a surface 
charge amount \9* times the membrane area 4ira 2 . To set the initial conditions 
for the Poisson model, we need to take into account this contribution. We 
spread this surface charge contribution uniformly throughout the intracellular 
and extracellular spaces. 

The problem with this initialization is that the excess charge should not be 
uniformly distributed but should be distributed so that the concentration profile 
shows a space charge layer near the membrane. Since we do not know the exact 
details of this concentration profile a priori, we let the Poisson system relax 
between — T r < ry < to a state where the bulk is approximately electroneutral 
and the excess charge accumulates near the membrane. During this period, 
we set the membrane current equal to zero and the dimensionless diffusion 
coefficients to be equal to Di = 1. We let T r = 10/3, 10 times the charge 
relaxation time. 

At time t = we turn on a current of constant strength a carried by ionic 
species i = 4 flowing from the extracellular space (r > a) into the intracellu- 
lar space (r < a). We let our simulations last until ry = T e = 2— , which 
is approximately the time it takes to depolarize the dimensionless membrane 
potential from —1 to 1. We place N r = 200 computational voxels in both the 
extracellular and intracellular regions (a total of 2N r — 400 voxels) , and we take 
the time step Ary = ^. Using a larger time step led to numerical instabilities 
with the Poisson model. A snapshot from a sample run of this simulation is 
shown in Figure [3] 

Raw data produced by the electroneutral model do not capture the ionic 
concentration or electrostatic potential profiles in the Debye layer. But it is 
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x=T /4 



x=3T /4 




Figure 3: Snapshots of simulation when (/?, a) = (10~ 35 , 10 -15 ). Three curves, 
the Poisson computation, the raw data and modified data from the electroneu- 
tral models arc plotted. The three curves arc virtually indistinguishable. 



28 



possible to produce an approximate profile in the Debye layer based on the 
asymptotic calculations we performed. We can see from (|73| and ([74]) that the 
Debye layer has the effect of adding a correction term to C, and $ that decays 
exponentially with distance from the membrane. The decay length and the 
magnitude of the correction term can be approximated by the values of $, d 
evaluated at the membrane and A^. For <£>, we modify the raw data from the 
clcctroncutral model as follows: 

^modified = $ + ( 153 ) 

5$ = _^exp(-r+^^) ifr>a (154) 
= J^ exp (- r -t_?t) iir<a (155) 

r± = 





yfz?Ci(r = a±) (156) 



where the double signs correspond in the last line. For the ionic concentrations 

^modified = (J, + SC . ( 15? ) 

A < (a+)g$ m r+ / r+ \r-a\ \ 

oCi = exp —1 T — - — it r > a (158) 

Zi V P J 

= exp —1 — - — if r < a (159) 

Zl V P ) 

We note that 8Ci and (5$ are expressed entirely in terms of raw data computed 
with the electroneutral model. When comparing the electroneutral model with 
the Poisson model, we shall use the above modified profile. 

In order to quantify the difference between the electroneutral and Poisson 
calculations, we introduce the following norm on the computational domain. 
Suppose the quantity u is defined at each voxel indexed by k. We define the 
discrete p-norm as: 



, El=l\v k \ , 



1 < p < oo (160) 



||w||^oo = max | life | (161) 

k 

where Uk is the value of u at the fc-th voxel and 14 is the volume of the fc-th 
voxel. In defining the LP norm in (|160|) . we have divided by a normalizing factor 
so that ||u|| iP gives an average measure of the "L p deviation". In particular, 
linip^oo ||u|| iP = ||w|| i0 o. For ionic concentrations Ci, we use the relative error: 



I /^electroneutral /^Poisson I 
I *~*4 — ^ i 

I ^Poisson | 



= — u7=^ ^ (162) 
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LP 


M p (d) 


M P {C 2 ) 


M P (C 3 ) 


M P (C 4 ) 


M P m 






3.79 x 10~ b 


6.98 x 10~ b 


6.56 x 10~ b 


1.55 x 10~ 4 


2.73 x 10~ b 


L' 2 


4.67 x lCT b 


8.43 x lQ- b 


7.71 x 10~ b 


1.74 x 10~ 4 


5.95 x 10~ b 




2.87 x 10~ 4 


2.35 x 10~ 4 


2.86 x 10~ 4 


5.85 x 10~ 4 


1.40 x 10~ 4 


(^2,Qt2) 


L v 


7.20 x 10~ b 


2.76 x 10~ 5 


3.25 x 10~ 5 


2.99 x 10~ 6 


4.70 x 10~ 5 


L' 2 


8.46 x 10~ b 


4.40 x 10~ b 


5.47 x 10~ b 


5.15 x 10~ b 


1.03 x 10~ 4 


L°° 


2.32 x 10~ 4 


1.57 x 10~ 4 


1.89 x 10" 4 


2.06 x 10~ 4 


1.71 x 10~ 4 


(^3,Q!3) 


L y 


1.74 x 10~ b 


1.79 x 10~ b 


4.25 x 10-'' 


3.31 x 10-' 


1.62 x 10~ b 


L' 2 


2.58 x 10"° 


2.98 x 10~ b 


1.97 x 10~ b 


1.71 x lO" 


4.16 x 10~ b 


L°° 


1.19 x 10~ b 


8.59 x lQ- b 


7.53 x 10~ b 


1.57 x 10~ 4 


6.95 x 10~ b 



Table 1: M. p values for spherical geometry for three computational experiments 
with different values of (3 and a. Here, (/3i,«i) = (10~ 3 ,10 -2 ), {(32, 0(2) = 

(10-3-5^0-1.5^ ^ 3ja3 ) = (1Q-4, 10-1). ' 



This is a more stringent criteria than using the absolute error (without the 
denominator in the above) especially for C4 whose initial concentration is very 
small. 

For the electrostatic potential $, it does not make sense to use the relative 
error since an arbitrary constant constant may be added to We thus, measure 
the error in $ as: 

Ep(§) = min ||$<^troncutral _ ^Poisson + ^JJ (163) 

Note that it is reasonable to consider the absolute error in <!>, since $ is di- 
mensionless, and its typical magnitude is 1. Though £ p ($) may in general be 
difficult to compute in closed form, this is possible when p = 1, 2,oo, values of 
p for which we shall compute £ P (&) in the following. 



In table |T]), we list the A4 p (u), u = $ or d where: 

MJu) = max SJu) (164) 

F 0<r v <T e F 

We see that for all parameter ranges tested here, the error falls within order 
10~ 4 . This translates to a 0.01% error in C, and an error of about 0.025mV 
in the dimensional electrostatic potential </>. In cases (f3,a) = (10~ 3 ,10~ 2 ) 
or (10~ 3 - 5 , 10~ 15 ), a is comparable in magnitude to -J]3~. The degree of cor- 
respondence exhibited for these two cases is remarkable since the asymptotic 
calculations were performed under the assumption that a = 0(1) with respect 
to a//3- It is notable that the relative error is order 10 -4 even for C4 which has 
a vanishing small concentration. This tells us that we may include ions of very 
small concentration into our model framework, which is significant if we are to 
include calcium dynamics pQ. 
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3D,Cj, 4> 


3D,Cj, 4> 




3D,0 


1D,0 






Electro- 










Cable 


Poisson 




neutral 




* 3D-Cable 






inner layer / 


int. layer 


outer layer/ 





Is thin domain asymptotics when a small 

matched asymptotics when (3 small 

Figure 4: Hierarchy of Electrophysiology Models 

We see that Moo{Ci) is significantly larger than jMi(Cj) or ./^(Ci)- Despite 
the modification we performed on the raw data for the clcctroncutral model, the 
deviation between the electroneutral and Poisson models are still concentrated 
at the Debye layer. Since this layer is very small in volume, the L 1 and L 2 
errors are not significantly affected. 

10 Conclusion 

The Poisson model, a candidate model for three dimensional cellular electrical 
activity, is computationally difficult to deal with, because of the presence of 
the Debye layer which develops at membrane interfaces. We introduced the 
electroneutral model as a computationally amenable and biophysically natural 
model of cellular electrical activity. We use asymptotic calculations to demon- 
strate the validity of the electroneutral model. The matched asymptotic cal- 
culations required the introduction of two boundary layers at the membrane, 
the inner Debye layer and the intermediate fast diffusion layer. We show that 
the electroneutral model gives an approximation to the Poisson model in the 
intermediate and outer layers as the small parameter (3, the ratio between the 
Debye length and the electrotonic length, becomes small. We demonstrated 
computationally that the electroneutral model gives an excellent approximation 
to the Poisson model. 

We have also succeeded in systematically deriving the standard cable model 
from the Poisson model or the electroneutral model. The above derivation can 
be viewed as a significant step toward a full study of the validity of the ca- 
ble model, an issue of fundamental importance to computational neuroscience 
|24j . In the course of this derivation, we have seen that there are models of in- 
termediate complexity in between the Poisson or electroneutral model and the 
cable model (Figure [4]) . The 3D-cable model and the simplified 3D-cable model 
describe the dynamics of the electrostatic potential in a three dimensional set- 
ting, but ignores the dynamics of ionic concentrations. We believe that each of 
these models will be suitable in certain situations, the Poisson or electroneutral 
models being the most detailed. 

When l m ~ vT? matching at the inner-intermediate layer interface resulted 
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in an additional surface drift-diffusion term along the membrane (cf. 155)) . It 
would be interesting to incorporate this into the electroneutral model and see 
whether this term leads to a significant difference in the behavior of the model. 
When l m ~ \f^ 1 matching between the intermediate and outer layers is probably 
challenging, since ionic fluxes parallel to the membrane will be comparable in 
magnitude to fluxes perpendicular to the membrane. The intermediate layer will 
lose its one-dimensional structure. We believe that the electroneutral model 
correctly captures the dynamics of ionic concentrations in the slow diffusion 
time scale (time scale T ). This claim is supported in part by the fact that the 
conservation relation, equation (|4ip is satisfied. We plan to investigate these 
points in future work. 



11 Appendix 

The calculations presented below are identical to the one that appears in [15], 
except for notational differences and some additions. We would like to solve ([55]) , 
([69]) under the boundary conditions ([70]) - ([72)1 . Since this is a one dimensional 
boundary value problem, we shall think of $ h0 and Cf as functions only of £j 
and do not explicitly write their dependence on ^\ or £3. 
Equation ([68]) can be integrated easily to obtain 



Cf (£ b ) = Cf (00) exp (-z 4 ($ h0 (^) - * b0 (oc))) . (165) 
This equation can be substituted into (|69[) to yield: 



<9 2 $ b0 ( 



Po 



J2 *Cf (00) exp (-*(* M (tf) - <P b0 (oo))) . (166) 



Here we use an approximation to linearize the above Poisson-Boltzmann equa- 
tion. We suppose 

|z i (* 60 (£ 1 6 )-$ W) (oo))| «1- (167) 

This can be justified if 9* is small, as was shown in [15]. The smallness of 
8* states that the amount of charge that may accumulate at the membrane is 
small. The smallness of this charge accumulation guarantees that the deviation 
of $ in the inner layer from its value in the intermediate layer is small. Given 
that (|167p is a valid assumption, we linearize (|166p to find: 

Cf (g) = Cf (00) (1 - *($ M (g) - $ b0 (oo))) (168) 
(* b0 (?i b )-* b0 (oo)) = r 2 ($ b0 (ei)-* b °(oo)) (169) 

N N 

r 2 = 5>?cf M = ^cfto), r>o. (170) 

»=X i=l 
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Here Cf°(0) is shorthand for Cf(tf = 0). To derive p^g|) and (US]), we have 
used 

JV N 

p + J2 ZiCf (oo) = Po + E = ( 171 ) 

i=l i=l 

which follows as a consequence of (|56[) and the matching condition (|71|) . Solving 
fn&ll with jMD and J75J), 

$ bo (^) = $ o0 (0)-^exp(-re l 6 ) (172) 

where $ a0 (0) is shorthand for $ o0 (^ = 0) and a is equal to 

a = <r($ bo (0)-$ (i) < bo (0)). (173) 
Hence, according to (|168|) and the matching condition (|71|) . 

C b0 (^) = Cf (0) (l + ^ expC-r^ 6 )) (174) 

We note that a is the total excess charge found in the inner layer, as can be 
seen as follows. The excess charge contributed by the i-th species of ion <Xj can 
be computed using expression (|174[) as: 

roc 2f<a0(n\ 

*i = J Zi(C?{&) Cf (0))c% b = p2 " = A^cr (175) 

From ([75]) . we see that Xi is given by: 

- _ zfCf(0) _ zfCfjO) 

We immediately conclude that A, = 1. The total excess charge is given 

by summing di in i. 

iv { N \ 

i=i \i=i / 

The factor A^ thus represents the fraction of excess charge contributed by the 
i-th species of ion. 

We now have the solutions C b0 and <& b0 except that a is expressed in terms 
of <& b0 . We shall now express a in terms of Cf °(0) and $ a0 (0). First, we observe 
by substituting £j = in (TTT2")) that 

a = r($ aO (0)-$ bo (0)). (178) 

We next rewrite $ a0 (0) = $( fc )- a0 (0), $ bo (0) = $( fc >- bo (0), 5 = ff( fc \r = r<» and 
consider (|173p and (|178[) as well as the corresponding expressions on the other 
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side of the membrane(the fi(" side). 

crW = 0*($( fc )> 6o (O) - $W- bo (0)) (179) 

j(fc) = r ( fe )($( fe )' o0 (0) - $W- b0 (0)) (180) 

erW = 0*($W' 6O (O) - $( fe )> bo (0)) (181) 

erW = rW(*W'°°(0) - $W« M (0)) (182) 
After some algebra, we find, 

a« = -er« = fl($W>°°(0) - $«' o0 (0)) (183) 
1 _ 1 1 1 

The meaning of relation (|184[) becomes clear once this is written in dimensional 
terms: 111 1 

c^ = ^ + 7^) + Wr (185) 

where T = jrd- This relation states that the effective membrane capacitance C m 
can be computed as the intrinsic membrane capacitance C*^ and the capacitance 
of the space charge layers eyW , cy(0 i n series. We note that in (|184[) . 8* is 
small in magnitude whereas r( fe ) and are order 1. Therefore, 6 sa 0* , and 
C m ~ C * . 
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In |15j , Aj was used in place of A* in (|38p . This expression substituted into the 
boundary condition ((55|) yields: 

ZjFi . n W = ^ +Q ] j ( i 86 ) 

By following the same procedure as in Section [7J it can be easily seen that 
(|3"2")) - (J3"4"]) together with (|186|) has the desired approximation properties. Unfor- 
tunately, this system is ill-posed. 

We shall exhibit the ill-posed behavior in a simple situation. Assume we have 
two regions, one intracellular and one extracellular. Let there be no transmem- 
brane currents. Suppose that there are two positive ionic species with identical 
physical properties: the valence and diffusion coefficient are equal and scaled to 
1. Assume moreover that the positive ionic charges are counterbalanced com- 
pletely by a spatially uniform immobile charge of magnitude — 1. Equations 
l[32" l) -([3i j) and (T86]) become: 

= ^ 1 +/3Vx-F i (187) 

OTV 

F< = - (V X C, + CiV x cf>) (188) 

1 = d + C 2 (189) 



34 



Let n be the outward normal pointing from intracellular to extracellular and the 
membrane potential [$] = — $ c . The boundary conditions on the intracellular 
and extracellular sides of the membrane are respectively: 



On 



aP = (190) 



dry 
*(«>) 

df ] = -c^e^] (i9i) 



d °i ' _ r c „ -(e) 



dry 

We solve the above with the following initial condition: 

[$](X,0) = $ ^0 (192) 
C;(X,0) = C,o(X),C ll0 (X) + C 2 ,o(X) = 1 (193) 

We thus assume that the membrane potential is initially constant(= $o) through- 
out, whereas the ionic concentration may be nonuniform. From a physical 
standpoint, the system should relax to an equilibrium state in which the ionic 
concentration gradients have disappeared. 

We now show that this initial value problem is ill-posed. Summing equation 
(|187[) in i and using (|189[) one immediately concludes: 

A x $ = (194) 

To obtain boundary conditions for the above Laplace equation, we take the 
summation of both (|190[) and (|191[) in i to obtain: 

dry on On 

The equations (| 194[) . p95[) . p92j) together form an initial value problem for $ 
and this has a unique solution: <f> docs not change, and is constant within each 
spatial region. 

We now turn to Cj. From equation (|187p we obtain: 

8C 

— = /3A x C, C = d, C 2 = l-d (196) 

OTy 

where we used $ = const within each region. The boundary conditions are: 

„ dC' 1 dC' 1 

= " 1ST < 197) 

dry on on c 

where n c = — n is the normal pointing from the extracellular to intracellular 
space. The evolution equations for the concentrations completely decouple into 
two separate diffusion problems for which the boundary conditions have the 
form kQpjr + ^ = 0. When k is negative, this problem is ill-posed, as was 
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formally established recently in [35]. We see from (|197[) and (|198p that one of 
the diffusion problems is bound to be ill-posed unless $o = identically. Here, 
we shall illustrate this by way of a simple example. 

Consider the above in X = (X, Y) £ M. 2 and let the upper and lower half 
planes correspond to the intracellular and extracellular spaces respectively. We 
let $o = — lj an d seek solutions to (|196[) and (|197[) in the upper half plane 
subject to the condition that C decays to as Y — > oo. We obtain a family of 
solutions parametrized by I > 1: 

C^ry) = exp [±.r v - ± Y ) sin (|v^) (199) 

If the initial data contain any non-zero frequency component along the mem- 
brane, this component will grow exponentially, the exponent being roughly pro- 
portional to the wave number. Thus, the problem is ill-posed. 

This instability is most probably a generic feature of the equations not con- 
fined to the simple situation above. The instability is caused by the J^- term 

in the boundary conditions, which came from the term. In general, the 
boundary conditions are complicated functions of the ionic concentrations, but 
the leading order terms and ^ will dominate in stability considerations. 

Since the membrane potential [$] is multiplying the term, the diffusion 
problem is bound to be ill-posed at least on one side of the membrane. 

We now take a closer look at the above situation in an attempt to obtain a 
well-posed system of equations. Equation (|199j) tells us that the time constant 
associated with exponential growth in the ill-posed solution is at most /39 2 , since 
I > 1 . This time duration belongs to the charge relaxation regime (actually even 
faster, by a factor of 9 2 ). The spatial scale that appears in (|199|) is on the order of 
the Debye length or shorter. The instabilities that develop are thus inconsistent 
with our ansatz that the evolution of Ci and $ do not possess spatiotemporal 
scales associated with charge relaxation in the space charge layer. 

We would like to remove the explosive behavior caused by . We propose 
the following fix. Let A 2 ; be a quantity that evolves according to the following 
differential equation. 

Thus Xi tracks A^ with a time lag t\ = /3, the charge relaxation time. This 
has the effect of filtering out any temporal structure that exists on a time scale 
smaller than 0((3). Instead of Aj, we shall use A; in (|186p . We note that since the 
relaxation time constant (= (3) is taken equal for all ionic species, the important 
relation ^ i A, = 1 holds true as long as this relation is satisfied at the initial 
time (see equation (|3D]l). 

It is important to demonstrate that this replacement does not change the 
formal approximation properties of the original system of equations. We can 
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(204) 



find the discrepancy between A 2 : and A; as follows. We can solve (|200[) so that: 
K(t v ) = Ai(0) exp (-T V IP) + ^J q V Ai(s) exp (^p 1 ) da (201) 
If ry is order 1, expanding A^(s) around ry , one can easily see that: 

Ai(Ty) =Ai(0) exp (-7V//9) + Ai(7v) (1 - exp (-r v /p)) 

q\. (202) 

+ /3 7r ^(Tv)(i-exp(-r y //?)) + ••■ 

CT v 

We see that 

Ai(7v)=Ai(7v) + 0(/3) (203) 
as long as is 0(1). Likewise, 

9A 1 <9A 

g-r(rv) = ~ ^A,(0) exp (-rv/P) + q^( 1 ~ cx P (-tv/0)) 

+ f3^(rv)(l-cxp(-r v /P)) + --- 
from which we find that 

^ {Tv) = T^ {Tv) + m (205) 

as long as is O(l). It is also possible to show that the error is 0(P) when 

t v = O(P) provided A 4 (0) = A;(0). Since A, and follow A l and to 
order O(P), replacing (|186[) with (|200p will not alter the formal approximation 
properties of the ill-posed model. 

Now we perform the same half plane analysis for the model we just proposed 
as was done for the ill-posed system. We take t\ in (|200|) as a parameter for 
now, and see what values of t\ will remove the instability. The expression 
corresponding to ()199|) is: 

C;(X v ,Ty) = exp (It v - mY a ) sm(kX a ) (206) 

k 2 = m 2 ~l, m =^L (207) 
Txl + l 

Exponential growth corresponds to I > 0. We would therefore like to make sure 
that the following equation for I does not have a positive solution for any real 

M) 2 -'- k2 {m 

This is equivalent to showing that the left hand side of the above is non positive 
when I > 0. Note that I > implies: 

I- m% l<(^ l)l (209) 
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Therefore, t\ = (3 is more than adequate to make the above expression non 
positive, since 9 is a small number much less than 1. We thus see that for the 
above situation in which model the system (|3"2")) - (j3"^l) . (|186|) fails, the new model 
is stable. 

What we have done is to add a stabilizing term to an asymptotically correct 
but ill-posed system. The situation here is analogous to having a consistent 
but unstable numerical discretization for an evolution equation. In such cases, 
one often adds to the numerical scheme a stabilizing term (e.g. small diffusive 
correction) whose order is small so that it does not alter the consistency of the 
scheme pH] . 
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